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ABSTRACT 

Observations of long-lived mixed-phase Arctic boundary layer clouds on 7 May 1998 during the First In- 
ternational Satellite Cloud Climatology Project (ISCCP) Regional Experiment (FIRE)-Arctic Cloud Ex- 
periment ( ACE)/Surface Heat Budget of the Arctic Ocean (SHEBA) campaign provide a unique opportunity 
to test understanding of cloud ice formation. Under the microphysically simple conditions observed (ap- 
parently negligible ice aggregation, sublimation, and multiplication), the only expected source of new ice 
crystals is activation of heterogeneous ice nuclei (IN) and the only sink is sedimentation. Large-eddy simu- 
lations with size-resolved microphysics are initialized with IN number concentration N in measured above 
cloud top, but details of IN activation behavior are unknown. If activated rapidly (in deposition, condensation, 
or immersion modes), as commonly assumed, IN are depleted from the well-mixed boundary layer within 
minutes. Quasi-equilibrium ice number concentration JV,- is then limited to a small fraction of overlying N in 
that is determined by the cloud-top entrainment rate w e divided by the number-weighted ice fall speed at the 
surface ty. Because w c < 1 cm s _1 and vr > 10 cm s -1 , NJN m -c 1. Such conditions may be common for this 
cloud type, which has implications for modeling IN diagnostically, interpreting measurements, and quanti- 
fying sensitivity to increasing N in (when w e hf < 1, entrainment rate limitations serve to buffer cloud system 
response). To reproduce observed ice crystal size distributions and cloud radar reflectivities with rapidly 
consumed IN in this case, the measured above-cloud N m must be multiplied by approximately 30. However, 
results are sensitive to assumed ice crystal properties not constrained by measurements. In addition, simu- 
lations do not reproduce the pronounced mesoscale heterogeneity in radar reflectivity that is observed. 


1. Introduction 

Observations indicate that the Arctic has warmed at 
roughly twice the global average rate since the preindustrial 
period, and that trend is expected to continue during this 
century (Solomon et al. 2007). However, climate model 
predictions vary considerably, owing at least in part to 
the complexity of atmosphere-ice-ocean interactions and 
a scarcity of the data required to study them (Randall et al. 
1998; Sorteberg et al. 2007). Differences in climate model 
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representation of clouds have been targeted as a cause 
for spread in Arctic climate predictions (Inoue et al. 
2006; Gorodetskaya et al. 2008; Holland et al. 2008). 

It is therefore a research objective to generate micro- 
physically detailed, high-resolution simulations of the most 
relevant cloud types in order to understand the dominant 
processes and improve their necessarily simplified rep- 
resentation in climate models. Because of the many 
known gaps in our knowledge of cloud processes, com- 
prehensive field experiment case studies are required to 
evaluate simulation fidelity. Here we consider an ob- 
served case of low-level mixed-phase clouds, a common 
and persistent cloud type over Arctic sea ice during the 
spring and autumn transition seasons (Shupe et al. 2006), 
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when sea ice is changing most rapidly in a manner that 
may be associated with cloud processes (e.g., Zhang et al. 
1996; Dong et al. 2001). This cloud type also appears to be 
particularly poorly represented in climate models owing 
at least in part to a lack of understanding of the relevant 
microphysical processes (Prenni et al. 2007). 

Of leading importance for constraining detailed sim- 
ulations of mixed-phase boundary layer clouds are in 
situ measurements of water droplet and ice crystal size 
distribution, ice crystal habit, and ice nucleus (IN) num- 
ber concentration iV IN active under in-cloud conditions. 
Ancillary meteorological measurements are required to 
provide model initial and boundary conditions, and 
ground-based cloud radar measurements provide valu- 
able additional constraints on model performance (e.g., 
Fan et al. 2009; van Diedenhoven et al. 2009). To our 
knowledge, only three field experiments to date have 
provided all such measurements for single-layer cases of 
shallow mixed-phase cloud that are most suitable for 
basic modeling case studies: the 1998 First International 
Satellite Cloud Climatology Project (ISCCP) Regional 
Experiment-Arctic Cloud Experiment (FIRE-ACE)/ 
Surface Heat Budget in the Arctic (SHEBA) campaign 
(Curry et al. 2000), the 2004 Mixed-Phase Arctic Cloud 
Experiment (M-PACE; Verlinde et al. 2007), and the 2008 
Indirect and Semi-Direct Aerosol Campaign (ISDAC; 
McFarquhar et al. 2011). 

Perhaps the most extensively studied measurements 
to date were obtained on 10 October during M-PACE in 
a supercooled boundary layer cloud (mixed-phase layer 
circa — 9° to — 16°C) that formed over the ice-free Beaufort 
Sea under clean, cold-air outbreak conditions (McFarquhar 
et al. 2007). In a broad model intercomparison study based 
on these observations and organized in association with the 
Global Energy and Water Cycle Experiment (GEWEX) 
Cloud System Study (GCSS) program, it was found that 
even high-resolution models with relatively sophisticated 
microphysics, when initialized and forced identically, 
produced widely differing results (Klein et al. 2009). 
Other studies of the case identified a controlling role for 
activated IN concentration in determining cloud prop- 
erties through the regulation of heterogeneous ice for- 
mation (Fridlind et al. 2007; Prenni et al. 2007; Morrison 
et al. 2008; Fan et al. 2009; Solomon et al. 2009; Avramov 
and Harrington 2010), consistent with analyses of ear- 
lier observed cases (e.g., Pinto 1998; Jiang et al. 2000). 
However, M-PACE modeling studies also reported large 
differences in the sensitivity of cloud properties to above- 
cloud N in , likely caused at least partly by differences in 
assumed ice crystal properties (Avramov and Harrington 
2010); as an aside we note that this sensitivity to above- 
cloud A in assumes Arctic boundary layer IN sources to 
be negligible (e.g., Pinto 1998; Harrington and Olsson 


2001). Finally, a subset of modeling studies concluded 
that A in measured above cloud were insufficient to ex- 
plain ice crystal number concentrations measured in the 
boundary layer (e.g., Fridlind et al. 2007; Morrison et al. 
2008; Fan et al. 2009), although it remains unknown 
whether large uncertainties assigned to observed ice 
crystal number concentrations (e.g., factor of 5; Fridlind 
et al. 2007) were adequate to account for errors associ- 
ated with ice crystal shattering on instrument probes (e.g., 
Korolev and Isaac 2005; Korolev et al. 2011). 

Modeling studies based on data gathered during the 
FIRE-ACE/SHEBA campaign (hereafter referred to as 
SHEBA) have also prominently identified the mecha- 
nisms of ice formation in mixed-phase boundary layer 
clouds as a leading source of uncertainty in model results 
(Girard and Curry 2001; Lohmann et al. 2001; Morrison 
and Pinto 2005; Morrison et al. 2005; Morrison and Pinto 
2006; Yuan et al. 2006; Sandvik et al. 2007; de Boer et al. 
2009). During SHEBA, observations of supercooled 
boundary layer clouds that formed under polluted con- 
ditions over sea ice on 7 May 1998 (mixed-phase layer 
circa —18° to — 20°C, droplet number concentration 
N d 200 cm~ 3 , and iV IN 2 L _1 ) provide a climato- 
logically important contrast to the 10 October M-PACE 
case of clean conditions over open ocean (N d «= 40 cm~ 3 
and X IN 0.2 L 1 ). The 7 May case has therefore been 
used as the basis for a follow-on model intercomparison 
study coordinated through GCSS (Morrison et al. 2011). 
From the standpoint of microphysical processes, the 
7 May GCSS SHEBA case is uniquely simple owing to 
high N d and relatively sparse concentrations of unrimed, 
nondendritic ice crystals, as discussed further below. 
This distinguishes it from the 10 October M-PACE case, 
with active drizzle and riming, and from the more re- 
cently observed 8 April ISDAC case, with active ag- 
gregation of dendrites (Avramov et al. 2011). 

Here we develop an adjusted version of the 12-h 
GCSS SHEBA case (see appendix) in order to better 
represent the last 2 h, when in situ ice particle size dis- 
tribution measurements were made. We use a large- 
eddy simulation code with size-resolved microphysics 
to simulate the coupling of dynamical and mixed-phase 
microphysical processes. Our principal objective is to 
determine whether the mean /V IN observed above cloud 
is adequate to explain mean observed boundary layer ice 
properties in simulations that are consistent with all 
other available observations. Since in situ measurements 
of ice crystal total number concentration were unreliable 
at sizes smaller than a poorly characterized threshold 
(e.g., Korolev et al. 2011), we compare simulations with 
1) in situ measurements of the size distribution of ice with 
maximum dimension larger than 200 /rm and 2) ground- 
based remote sensing measurements of cloud radar 
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reflectivity and mean Doppler velocity. Below we first 
describe the observations (section 2) and the model (sec- 
tion 3). We present a range of simulations using several 
approaches to represent IN and compare results with 
observations (section 4). Conclusions and implications are 
then summarized (section 5). 

2. Observations 

The 7 May 1998 flight of the National Center for At- 
mospheric Research (NCAR) C-130 aircraft was, to our 
knowledge, the only flight over SHEBA surface in- 
struments that took place in a long-lived (>12 h) mixed- 
phase boundary layer cloud deck without the overlying 
cloud layers that were commonly present (Wylie 2001). 
From the aircraft, we use measurements from a cloud 
particle imager (CPI), forward scatter spectrometer probe 
(FSSP-100), and two-dimensional cloud (2D-C) optical 
array probe (Lawson et al. 2001; Zuidema et al. 2005; 
Lawson and Zuidema 2009). We adopt the analysis of 
aerosol and IN data prepared for the GCSS SHEBA case 
(Morrison et al. 2011), which was based on aircraft mea- 
surements from a condensation nucleus counter (Yum 
and Hudson 2001) and counterflow diffusion chamber 
(CFDC; Rogers et al. 2001; Prenni et al. 2009). We use 6-h 
soundings and hourly surface measurements compiled for 
use by modelers (Persson et al. 2002; Beesley et al. 2000), 
which include liquid water path derived from microwave 
radiometer measurements (Liljegren 2000). We derive 
large-scale forcings from the National Centers for Envi- 
ronmental Prediction (NCEP)/NCAR 40-yr reanalysis 
project (Kalnay et al. 1996). We use radar reflectivity and 
mean Doppler velocity measurements from a K a -band 
millimeter-wavelength cloud radar (MMCR; Shupe et al. 
2001; Intrieri et al. 2002; Shupe et al. 2006). 

In brief overview, on 7 May 1998 the SHEBA ice 
station was located at roughly 75°N, 165°W beneath 
a widespread boundary layer cloud deck (Fig. 1) advect- 
ing to the northeast at about 5 m s 1 . During the 1200- 
2400 UTC period of the GCSS SHEBA case, MMCR 
measurements indicate cloud top decreasing from roughly 
600 to 400 m (Morrison et al. 2011, their Fig. 2). Aircraft 
measurements were limited to the last 2 h of this time pe- 
riod, 2200-2400 UTC. Several passes were made through 
the cloud layer (cloud droplets present), and two longer 
legs sampled ice properties beneath cloud base (Fig. 2). 

Although the best available aircraft altitude data in- 
dicate unphysically low elevations during a short period 
of the near-surface leg, here we use the altitude data 
only to separate particle size distribution measurements 
into in-cloud and below-cloud categories. At reported 
altitudes of 310-430 m, all FSSP concentrations indicate 
highly peaked droplet size distributions that were by 



Fig. 1. Advanced Very High Resolution Radiometer (AVHRR) 
channel 4 (10.5-12 /im) infrared satellite image at 2219 UTC 7 May 
1998. Figure reproduced from experiment Web site. 


contrast absent below 280 m (Fig. 3), thus indicating 
a cloud base range of 280-310 m that is reasonably 
consistent with ground-based lidar measurements (not 
shown). We use FSSP measurements only during these 
in-cloud time periods and only for diameters less than 
20 jum (Fig. 3a). We use ice measurements only below 
cloud base, noting that ice properties typically vary little 
with elevation in this cloud type (e.g., McFarquhar et al. 
2007, 2011). We also use 2D-C data only at maximum 
dimensions greater than 200 gm (Fig. 3b), where shat- 
tering effects on number concentrations could be less than 
about 20% at the small characteristic ice particle sizes 
observed here (Field et al. 2006). However, owing to the 
high degree of uncertainty associated with all optical array 
probe measurements (Korolev and Isaac 2005; Korolev 
et al. 2011), we perform an integrated analysis of the in 
situ and remote sensing measurements. 



Time (UTC) 


Fig. 2. Reported near-surface elevation of the C-130 aircraft 
during 2200-2400 UTC 7 May 1998 (solid line). Based on FSSP 
measurements, flight times are identified as in cloud (310-430 m 
bounded by dashed lines) and below cloud base (280 m indicated by 
dotted line). 
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Fig. 3. Observed hydrometeor size distributions measured with 
FSSP and 2D-C probes during 2229-2359 UTC 7 May 1998 (a) at 
310-430 m, (b) below 280 m, and (c) at all altitudes. Mean distri- 
butions in (a) and (b) shown in black solid lines for limited size 
ranges are reproduced in (c) for comparison with mean values at all 
elevations shown in black dashed line for all sizes. 


3. Model description 

a. Dynamics 

We use the Distributed Hydrodynamic Aerosol and 
Radiative Modeling Application (DHARMA) code, 
which treats dynamics using a large-eddy simulation (LES) 
model (Stevens et al. 2002). We use a horizontal domain 
that is 3.2 km on a side, roughly 7 times the boundary layer 
depth. A vertical extent of 1 km allows boundary layer 


depth evolution that is not affected by damping of gravity 
waves above 800 m through relaxation of potential tem- 
perature and horizontal winds toward their time-varying 
horizontal averages with a time scale of 100 s, applied at 
full strength at the domain top and decreasing as sine- 
squared to zero at 800 m. Grid spacing is uniform hor- 
izontally (50 m) and vertically (10 m). A dynamic 
Smagorinsky model (Kirkpatrick et al. 2006) is used to 
compute subgrid-scale mixing. Surface turbulent fluxes 
are computed from Monin-Obukhov similarity theory 
using the dimensionless profiles of Businger et al. (1971) 
with a turbulent Prandtl number of unity and a von 
Karman constant of 0.41. Skin water vapor is assumed 
saturated with respect to ice at the skin temperature. A 
surface roughness of 0.4 mm is assumed for momentum, 
water vapor, and heat (cf. Brunke et al. 2006). Hori- 
zontal winds are nudged toward their initial profiles with 
a 1-h time scale. The domain is translated with mean 
cloud-layer winds (1.8 m s _1 westerly and 4.3 m s _1 
southerly) to minimize errors associated with advection 
and allow vertical wind speed to dictate the maximum 
advective Courant number. A 5-s dynamical time step 
is taken unless the Courant number exceeds 0.8 (the 
strongest vertical wind speeds increase the total number 
of time steps by about 5% in a typical simulation here). 
Halving vertical and horizontal grid spacing to 5 and 
25 m, respectively, and halving the maximum time step 
to 2.5 s decreases cloud-top entrainment and increases 
liquid water path by about 10%, suggesting that the 
baseline resolution allows for reasonable representation 
of boundary layer dynamics when using DHARMA for 
this case. 

The specified profile of horizontally uniform large- 
scale subsidence is treated separately from the resolved 
vertical winds and only appears through a source term 
for each prognostic variable <f>, computed through first- 
order upwind advection as — w LS = — 5</>/5z, where z is 
altitude (cf. Wyant et al. 1997; Ackerman et al. 2009). 
We note that the cloud-top entrainment rate is com- 
puted throughout as the sum of the subsidence rate 
at the mean height of the boundary layer top (which is 
the same as cloud top here) plus the rate of change of 
mean boundary layer depth (cf. Faloona et al. 2005, their 
Eq. 2). 

b. Microphysics 

We use size-resolved microphysics based on the 
Community Aerosol-Radiation-Microphysics Applica- 
tion (CARMA) code. The microphysical formulations for 
warm and cold clouds are described by Ackerman et al. 
(1995) and Jensen et al. (1998), respectively. An earlier 
version of the mixed-phase formulation is described by 
Fridlind et al. (2007), and modifications since that study 
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are described below. The linkages between microphysics 
and dynamics are described by McFarlane et al. (2002, 
their appendix B). 

We use 32 mass-doubling bins each for droplets and 
ice, where the mass of the smallest bin in each grid is that 
of a droplet with diameter 2 pm. The mass of the largest 
bin is set by the requirement that it contain negligible ice 
under simulated conditions. Time substepping is em- 
ployed with a minimum step of 0.2 s to locally resolve 
fast microphysical processes such as droplet activation 
and condensational growth. Aerosols are initialized as 
specified in the GCSS SHEBA case and treated di- 
agnostically (Clark 1974) to avoid the need for 1) aerosol 
source terms, which are unknown, and 2) core second 
moments to restore aerosol size dispersion upon droplet 
evaporation (Ackerman et al. 1995), which are compu- 
tationally expensive. 

We treat ice in each size bin using the approach de- 
veloped by Bohm (1989, 1992a,b,c, 1994, 1999, 2004), 
which provides an integrated treatment of fall speeds 
and collision efficiencies for ice and liquid particles based 
on four properties of each participating particle type: 
mass, maximum dimension, projected area, and aspect 
ratio. We use size-dependent coalescence efficiencies for 
water droplets (Beard and Ochs 1984), a coalescence ef- 
ficiency of unity for liquid-ice collisions and 0.1 for ice-ice 
collisions of nondendritic crystals under dry-growth 
conditions (e.g., Mitchell 1988; Wang and Chang 1993). 
Results are negligibly impacted by increasing the ice-ice 
collision efficiency to 0.25 or 0.3 (e.g., Mitchell 1988; 
Girard and Blanchet 2001). We neglect any turbulence 
effects on the gravitational collection process, which are 
likely to be minimal under the relatively weak dynami- 
cal conditions of this case. 

To approximate the impact of ice habit on vapor de- 
position and sublimation rates, capacitance is calculated 
for oblate spheroids [Pruppacher and Klett 1997, their 
Eq. (13-78)], where the aspect ratio is taken as the ratio 
of minor to major axis and major axis is maximum di- 
mension. We consider the impact of reduced capacitance 
in a sensitivity test below: 1) the ratio of capacitance to 
maximum dimension is specified to be 0.35 for all particle 
sizes [Westbrook et al. 2008, plate aspect ratio of 0.1 in 
their Eq. (3)], and 2) particles of 120-240-/xm maximum 
dimension are assumed to comprise a linearly increasing 
fraction of aggregates with a ratio of capacitance to max- 
imum dimension reduced to 0.25 (Westbrook et al. 2008). 

c. Ice properties 

Ice is commonly represented in microphysics schemes 
by a fixed number of types such as plates or dendrites with 
predetermined properties that are not varied on a case- 
specific basis (e.g., Lynn et al. 2005). Less commonly, ice 


properties may be dynamically predicted (e.g., Morrison 
and Grabowski 2008; Hashino and Tripoli 2007), allowing 
case-specific properties to emerge in simulations. How- 
ever, since ice properties vary significantly even within 
basic habit classes, and evaluating the prediction of ice 
crystal habit is not an objective of this study, here we 
choose case-specific model settings to represent the ob- 
served ice properties. Since the ice properties needed by 
this model are not directly measured (viz., maximum di- 
ameter, projected area, and aspect ratio as a function 
of ice particle mass), the remainder of this section pro- 
vides an analysis of observations to derive ice properties 
consistent with the available CPI, 2D-C, and MMCR 
measurements. 

Manual examination of the available CPI images in- 
dicates an array of crystal shapes (e.g., Fig. 4) that is 
consistent with past observations at —16° to — 20°C (e.g., 
Magono and Lee 1966; Korolev et al. 1999) and labo- 
ratory observations of ice grown at those temperatures 
under conditions of 10%-20% ice supersaturation (Bailey 
and Hallett 2002; Bacon et al. 2003; Bailey and Hallett 
2004). A minority are relatively pristine plates with some 
degree of transparency (habit class Pla). A few plates 
have sectorlike branches, consistent with the warm end of 
the boundary layer temperature range (Plb; Magono and 
Lee 1966) or are nonisometric (e.g., Magono and Lee 
1966; Bacon et al. 2003). Most ice crystals appear poly- 
crystalline, including plates with spatial sectors (P5a) and 
radiating assemblages of plates (P6a). Some are small 
assemblages of minute plates or irregular germs (G5 and 
G6; Magono and Lee 1966; Bailey and Hallett 2004). 
Many larger crystals are what Bailey and Hallett (2004, 
p. 521) refer to as “jumbled arrangements of poorly 
formed but faceted plates or polyhedra of nonhexagonal 
shape” that at high supersaturation appear as “spatially 
extended forms.” 

Historically it has been common to express the re- 
lationship of particle maximum dimension D to mass 
m through power laws of the form m = aD h . We use 
several mass-dimensional relationships to span our ice 
mass grid piecewise. For instance, when ice crystals with 
D > 120 pm are assumed to be radiating assemblages of 
plates, ice with D <5 pm is treated as spherical, and ice in 
the transition size range is represented based on a power 
law transition between the mass of a sphere with D = 5 ftm 
and the mass of a radiating plate with D = 120 pm 
(Table 1). To choose a baseline mass-dimensional rela- 
tionship for the largest particles, we also considered two 
other habit choices based on ice shapes seen in CPI 
images: pristine hexagonal plates and aggregates that 
include plates (see Table 1). As an observation-based 
test of the validity of each candidate relation, we com- 
pare MMCR measurements with 35-GHz reflectivities 
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Fig. 4. Observed CPI images of ice crystals at 2301 UTC 7 May 1998 (aircraft below cloud at 
height of —200 m; cf. Fig. 2), representative of ice sampled over 2200-2400 UTC. Length scale 
shown at top. 


calculated from the individual in situ ice size distribu- 
tions below cloud shown in Fig. 3 (all particle sizes ini- 
tially included). Following the method described by van 
Diedenhoven et al. (2009), measured ice particle size 
distributions were averaged over 30-s time periods, each 
mass-dimensional relation assumed for particles with 
D > 150 /ttm (smaller particles assumed spherical), and 
reflectivities calculated using the QuickBeam package 
(Flaynes et al. 2007). Comparison with the available 
below-cloud MMCR reflectivity measurements in the 
same time ranges (as a proxy for the same locations) 
indicates that assuming large particles are radiating 


assemblages of plates results in calculated reflectivities 
that agree best with observations (Fig. 5). Particles of 
D < 200 gm do not contribute significantly (Fig. 5d). We 
do not consider this test particularly robust owing to the 
variability of MMCR reflectivity with time and the 
relatively sparse aircraft sampling. Nonetheless, ow- 
ing to a lack of other constraints, we assume radiating 
plates to represent large ice, adopt the ad hoc area- 
dimensional relationship proposed by Mitchell (1996), 
and perform sensitivity tests below. 

The relationship of aspect ratio to maximum di- 
mension has not received as much attention as mass and 
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TABLE 1. Mass- and area-dimensional power laws used in simulations and radar reflectivity calculations. 


Habit 

D (cm) 1 

a b 

b b 

c b 

d b 

Source 0 

Spheres 

0.0002-0.0005 

0.480 14 

3.00 

0.785 40 

2.00 

— 

Transitional 

0.0005-0.012 

0.023 06 

2.61 

0.175 96 

1.82 

— 

Radiating assemblages of plates 

>0.012 

0.002 40 

2.1 

0.228 50 

1.88 

LH74, MZP90, 

Aggregates of unrimed radiating 

>0.012 

0.002 94 

1.9 

0.228 50 

1.88 

M96, BL06 
LH74, M96 

assemblages of plates, side planes, 
bullets, and columns 
Plates with sectorlike branches 

0.001-0.016 

0.006 14 

2.42 

0.24 

1.85 

M96 


>0.016 

0.001 42 

2.02 

0.55 

1.97 

M96 

Hexagonal plates 

>0.012 

0.007 39 

2.45 

0.65 

2.00 

M96 


a Range of maximum crystal dimension D over which relationships are applied piecewise in simulations. Ranges shown for spheres and 
transitional properties (see section 3c) are those used when the largest ice crystals are radiating assemblages of plates. When the largest 
crystals are aggregates, the properties of spheres are applied over 0.0002-0.012 cm. When the largest crystals are sectored plates, the 
properties of spheres are applied over 0.0002-0.001 cm and the two consecutive relations shown are then applied. Hexagonal plates are 
used only in radar reflectivity calculations (cf. section 4d). 

b Values of a , b , c, and d in mass- and area-dimensional power laws m = aD h and A = cD d , where m is mass in grams, D is maximum 
dimension in centimeters, and A is projected area in centimeters squared. 
c LH74 = Locatelli and Hobbs (1974), MZP90 = Mitchell et al. (1990), M96 = Mitchell (1996), BL06 = Baker and Lawson (2006). 


area. Based on an analysis of aspect ratio using CPI 
images collected in Arctic clouds in the —15° to — 20°C 
range (Korolev and Isaac 2003), we assume that the as- 
pect ratio decreases linearly from 1.0 to 0.6 over a maxi- 
mum dimension range of 5-120 gm and remains constant 
at larger sizes. Aspect ratio primarily influences ice fall 
speed and capacitance (see sensitivity tests in section 4d). 

d. Ice formation 

Since all ice crystals are present in the boundary layer 
temperature range of —16° to — 20°C and no ice was 
observed to be seeding the cloud from above, we assume 
that all primary ice nucleation proceeds heterogeneously. 
We take two approaches to represent heterogeneous IN 
activation. First, we follow the simplified diagnostic ap- 
proach specified in the GCSS SHEBA case. The IN are 
activated as ice crystals if ice supersaturation exceeds 5% 
and are added to each grid cell such that the sum of ice 
crystals and IN never falls below the initial concentration 
of IN. A diagnostic approach was selected for the in- 
tercomparison based on results of the 10 October 
M-PACE model intercomparison, in which predicted 
ice crystal number concentrations ranged over five or- 
ders of magnitude, and it was therefore recommended 
that future studies constrain the treatment of ice nucle- 
ation and ice crystal number concentration (Klein et al. 
2009). However, since ice crystal number concentration is 
itself highly uncertain, as are virtually all details of IN 
activity, this simplified approach is based on the sugges- 
tion that ice concentrations are roughly equal to over- 
lying IN concentrations in this cloud type (Prenni et al. 
2007). A general consequence of the diagnostic approach 
is that any IN consumption is compensated by an 


unlimited source of IN replenishment (e.g., Harrington 
and Olsson 2001). 

We alternatively take a prognostic approach that ac- 
counts for IN sources, consumption, and transport 
(Fridlind et al. 2007). A spectrum of IN in each model 
grid cell is tracked in an array that ranges from least to 
most easily nucleated. Each array member contains IN 
that could be activated in any of the four commonly 
accepted modes: deposition, condensation, immersion, 
and contact (e.g., Pruppacher and Klett 1997). To cal- 
culate the rate of scavenging in the contact mode, all IN 
are assumed to be 0.5 gim in diameter, the mean effec- 
tive dimension observed during SHEBA (Rogers et al. 
2001). We assume that sublimated ice crystals yield IN 
that are preactivated (e.g., Roberts and Hallett 1968; 
Knopf and Koop 2006) and therefore in the array 
member that is easiest to nucleate, but alternatively 
assuming no IN regeneration from sublimated crystals 
changes results negligibly since the air is saturated with 
respect to ice nearly to the surface in this case. Generic 
IN activation properties are assumed (Fridlind et al. 
2007, their Table 1); they are not readily obtained from 
CFDC measurements because the instrument is not 
designed to distinguish between modes of activation and 
high spatial variability is commonly encountered during 
instrument scans over operating conditions in-flight (e.g., 
Rogers et al. 2001; Prenni et al. 2007). To represent ob- 
served conditions in this study, we initialize /V IN to 1.7 L 1 
based on the analysis of CFDC measurements conducted 
for the GCSS SHEBA case. Given the generic IN acti- 
vation properties (e.g., IN availability increases linearly 
in the condensation mode over the temperature range 
—8° to — 22°C), all 1.7 L are accessible under boundary 
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Fig. 5. Observed and calculated radar reflectivities. MMCR reflectivity observed between the surface and cloud 
base (180-280 m) during measurement of in situ size distributions over 2200-2400 UTC 7 May 1998 (shaded with 
median in dashed white line). Radar reflectivities calculated from size distributions observed in situ use varying mass- 
dimensional relations (solid lines with median in dashed black line; see section 3c and Table 1). 


layer conditions in two modes (deposition and contact). 
For sensitivity tests in which all IN operate in only one 
mode at a time, the condensation mode temperature range 
is limited to —8° to — 19°C and the immersion mode 
temperature range is limited to —10° to — 19°C; this 
guarantees that all IN can be activated in each mode 
independently under in-cloud conditions. 

No well-established ice multiplication processes ap- 
pear capable of significant secondary ice production 
under the observed conditions. Since liquid water con- 
tent is small and was found only at temperatures colder 
than — 18°C, Flallett-Mossop rime splintering is not ac- 
tive (Heymsfield and Mossop 1984). Shattering of drops 
larger than 50 gm in diameter is included as described by 
Fridlind et al. (2007), but the simulated number con- 
centration of such large drops is too small to be relevant. 
Simulated ice splinter production via ice-ice collisions is 
found to be insubstantial here when adopting an upper 
limit on the likely rate using the Vardiman (1978) 


parameterization as described in Fridlind et al. (2007), 
although this may not represent the maximum possible 
source because the unknown degree of ice crystal fall 
speed diversity is underestimated by choosing a single 
set of properties for each ice mass bin. Aside we note 
that such fall speed diversity was also neglected in our 
simulations of the 10 October M-PACE case (Fridlind 
et al. 2007), where observations indicated the coexistence 
of rimed and dendritic ice types likely more conducive 
to such multiplication (cf. Vardiman 1978; Yano and 
Phillips 2011). 

e. Radiative transfer 

Radiative transfer in 44 wavelength bins is computed 
independently for each column every 60 s using a two- 
stream model (Toon et al. 1989) in which the water vapor 
continuum absorption has been modified (Clough et al. 
1989). Particle scattering and absorption coefficients are 
computed from Lorenz-Mie theory (Toon and Ackerman 
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1981). Since longwave fluxes outside of the 4.5-62-/rm 
wavelength range are not rigorously included (their im- 
pact on simulations is negligible), we account for their 
contribution when comparing with measurements by 
adding 6.7 W m~ 2 (the average flux in that wavelength 
range under simulated conditions). For radiative transfer, 
ice is treated as spherical with diameter equal to maxi- 
mum dimension; this will be improved in future model 
development, but is sufficient for this case since ice has 
little impact on radiative fluxes (see also Zuidema et al. 
2005). 

4. Results 

a. Model setup 

We initialize model thermodynamic profiles, surface 
conditions, and top-of-model downwelling radiative fluxes, 
and apply large-scale tendencies over the 4-h simulation 
duration based on our adjustment of the GCSS SFIEBA 
case (see appendix). As in the baseline GCSS SFIEBA 
case, N m is initialized to 1 .7 L 1 and aerosol are initial- 
ized in two lognormal modes with geometric standard 
deviations of 2.04 and 2.5, geometric radii of 0.052 and 
1.3 pm, and number concentrations of 350 and 2 cm~ 3 , 
respectively. When ice is not present, these conditions 
lead to a cloud-topped boundary layer with steady liquid 
water path (LWP; see appendix). 

b. Diagnostic versus prognostic IN 

To introduce ice formation, we first use a diagnostic 
treatment of IN, which sustains the ice crystal concen- 
tration continuously at the initial Ain of 1.7 L 1 (see 
section 3d). This results in complete desiccation of the 
initial liquid water cloud within the 4-h simulation time 
(Fig. 6, solid lines). We next use a prognostic treatment 
of IN, which accounts for IN consumption (Fig. 6, dotted 
lines). After initial boundary layer IN are quickly con- 
sumed, the only (weak) source of new IN is then cloud- 
top entrainment, and LWP reaches a quasi-equilibrium 
state (defined throughout as sustaining an e-folding life- 
time of at least 10-h during hours 3^1; see appendix). 
Most IN are consumed instantly, and once boundary 
layer turbulence develops, the remainder are consumed 
within minutes. Since the boundary layer is saturated 
with respect to ice in this case, sublimation is negligible 
and activated IN are removed from the boundary layer 
when ice crystals sediment. 

If the simulation with prognostic IN is repeated with 
all IN available in only one mode at a time (see section 
3d), then consumption remains similarly efficient for all 
modes except contact, as discussed further below. The IN 
are activated at a much slower rate in the contact mode 



Fig. 6. Simulated domain-mean LWP, boundary layer depth 
H (defined by mean elevation where liquid water potential tem- 
perature is 258 K), droplet number concentration N,i (averaged 
over all grid cells with liquid water mixing ratio >10 -3 g kg -1 ), 
boundary layer (BL) IN and ice crystals iV,- averaged over H, 
domain-mean ice water path (IWP), domain-maximum variance of 
vertical wind speed W, and cloud-top entrainment rate w e (com- 
puted as dHIdt plus large-scale subsidence rate at H). Simulations 
listed in Table 2: diagnostic IN (solid lines), prognostic IN (dotted 
lines), contact IN only (short dashed lines), steady-state prognostic 
IN (dash-dotted lines), baseline (dash-triple-dotted lines), and 
IN X 30 (long dashed lines). Boundary layer IN and IV,- remain 
small at all times in some simulations. 

owing to inefficient scavenging of IN by droplets, and are 
scarcely consumed from the boundary layer within the 
4-h simulation time (Fig. 6, short dashed lines). Since it is 
not expected that IN active in the contact mode are in- 
active in other modes (e.g., Prenni et al. 2009), and we 
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TABLE 2. Simulations with size-resolved mixed-phase microphysics. 


Simulation 

Case 

specification 3 

IN scheme 

Initial 

above-cloud N m (L _1 ) 

Initial 

boundary layer N in (L~') b 

Ice crystal habit c 

Diagnostic IN 

A 

Diagnostic 

1.7 

1.7 

Radiating plates 

Prognostic IN 

A 

Prognostic 

1.7 

1.7 

Radiating plates 

Contact IN only 

A 

Prognostic 

1.7 

1.7 

Radiating plates 

Steady-state 

A 

Prognostic 

1.7 

0 

Radiating plates 

prognostic IN 
Baseline 

B 

Prognostic 

1.7 

0 

Radiating plates 

IN X 30 

B 

Prognostic 

51 

0 

Radiating plates 

Deposition IN only 

B 

Prognostic 

51 

0 

Radiating plates 

Condensation IN only 

B 

Prognostic 

51 

0 

Radiating plates 

Immersion IN only 

B 

Prognostic 

51 

0 

Radiating plates 

Decreased capacitance 

B 

Prognostic 

51 

0 

Radiating plates 

Aggregates 

B 

Prognostic 

51 

0 

Aggregates with 

Plates 

B 

Prognostic 

51 

0 

plates 

Sectored plates 

Modified diagnostic IN 

B 

Diagnostic 

0.29 

0.29 

Radiating plates 

GCSS submission 

G 

Diagnostic 

1.7 

1.7 

Radiating plates 


a G = original GCSS SHEBA case specification for 1200-2400 UTC (Morrison et al. 2011), A = adjusted case specification for 2000- 
2400 UTC (see appendix), B = baseline case specification with increased downwelling longwave radiation and moisture convergence 
(see section 4c). 

b Initializing boundary layer N m to zero reduces spinup associated with starting far from quasi-equilibrium (see section 4b). 
c Habit of ice crystals with largest maximum dimensions; smaller particles spherical or transitional (see section 3c and Table 1). 


have no evidence for an independent reservoir of 
contact IN, we assume that any contact IN can act in 
at least one other mode and are therefore activated 
rapidly. Based on the first published simulations that 
apply a prognostic approach to ice nucleation in mixed- 
phase boundary layer clouds, Harrington and Olsson 
(2001) also describe rapid IN depletion. Others have 
reported it in simulations of the 10 October M-PACE 
case (Fridlind et al. 2007; Fan et al. 2009; Avramov and 
Harrington 2010) and under SHEBA conditions (Morrison 
et al. 2005). 

Whereas all liquid water was consumed when treating 
IN diagnostically, desiccation is by contrast limited 
when treating IN prognostically, despite an initial burst 
of ice formation that does not persist when accounting 
for IN depletion. To eliminate the initial burst of ice 
formation and more quickly reach quasi-equilibrium 
ice water path (see also Fridlind et al. 2007), we next 
initialize IN in the boundary layer to zero, leaving only 
IN above the boundary layer at the background value 
of 1.7 L 1 (Fig. 6, dash-dotted lines). We refer to this 
as a steady-state initialization approach (since bound- 
ary layer IN concentration is initialized close to its 
very low quasi-equilibrium value), and use il in the re- 
maining simulations with prognostic IN (Table 2). 

c. IN insufficient to explain observed ice 

The simulated droplet number size distributions match 
in-cloud observations quite well in the simulation with 
steady-state prognostic IN, but the predicted number 


concentration of ice is too low by more than an order of 
magnitude at all sizes (Fig. 7a). Thus, N iN measured 
above cloud appears insufficient to explain observed ice 
crystal numbers. Using the QuickBeam package (Haynes 
et al. 2007) to calculate 35-GHz reflectivities and mean 
Doppler velocities below cloud from simulated ice crystal 
size distributions and vertical wind speeds at degraded 
vertical model resolution to match MMCR observations, 
as described by van Diedenhoven et al. (2009), we find 
that median simulated radar reflectivity is also greater 
than 10 dBZ lower than observed (Fig. 8a). This dis- 
crepancy is consistent with the model underestimation 
of ice number concentration over all observed sizes. Al- 
though the median of mean Doppler velocities is under- 
estimated by about 10 cm s 1 relative to the observed 
median of 50 cm s _1 (Fig. 8b), it is estimated that mea- 
sured Doppler velocities are biased high by about 
10 cm s 1 based on the shipborne radar tilt and boundary 
layer winds during 2000-2400 UTC. Thus agreement of 
the observed and simulated medians appears close, but 
the simulated distribution of mean Doppler velocities is 
broader than observed. 

The overly broad distribution of Doppler velocities 
suggests that simulated boundary layer dynamics may be 
too strong. Given the limitations of our modeling ap- 
proach and the constraints imposed by observed sur- 
face and sounding measurements, we are left with few 
relevant degrees of freedom. If the downwelling long- 
wave radiative flux specified at 1-km height, which is 
not directly constrained by observations in the GCSS 
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Fig. 7. Observed and simulated droplet and ice particle size distributions. Observed mean size distributions 
measured with FSSP and 2D-C probes during 2229-2359 UTC 7 May 1998 (solid black lines, as in Fig. 3) at reported 
aircraft altitudes of 310-430 m (drops) and below 280 m (ice) are compared with simulated size distributions of drops 
at 310-430 m and ice below 280 m (gray lines, black dashed line is mean). Simulations listed in Table 2: (a) steady- 
state prognostic IN, (b) baseline, (c) IN X 30, (d) aggregates, (e) plates, (f) decreased capacitance, (g) modified 
diagnostic IN, and (h) GCSS submission. Simulations are sampled at 12 h (GCSS submission, corresponding to 
2400 UTC) or 3 h (all others, corresponding to 2300 UTC). 


SHEBA case, is increased by 15 W m , then cloud- 
top radiative cooling and entrainment are reduced. 
If large-scale horizontal advective flux convergence 
of the water vapor mixing ratio q v is increased to a 
vertically uniform rate of 0.09 g kg -1 day -1 to maintain 


quasi-equilibrium LWP, simulations remain consistent 
with observations (see appendix). These changes result 
in the mean Doppler velocities agreeing better with 
measurements (Fig. 8d), with little associated impact 
on reflectivity (Fig. 8c versus Fig. 8a) or ice size 
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Fig. 8. Observed and simulated histograms of radar reflectivity and mean Doppler velocity below cloud base (180- 
280 m). Observed MMCR reflectivity and Doppler velocity during (a)-(l) 2200-2400 UTC and (m),(n) 1200- 
2400 UTC 7 May 1998 (shaded, dashed white line is median). Simulations listed in Table 2: (a),(b) steady-state 
prognostic IN, (c),(d) baseline, (e),(f) IN X 30, (g),(h) aggregates, (i),(j) plates, (k),(l) modified diagnostic IN, and 
(m),(n) GCSS submission. Simulations are randomly sampled (solid black line, dashed black line is median) over 
2-12 h (GCSS submission, corresponding to 1400-2400 UTC) or 3-4 h (all others, corresponding to 2300-2400 UTC). 
MMCR Doppler velocity is likely biased high by —10 cm s _1 based on radar tilt and boundary layer wind speeds 
during 2000-2400 UTC. 


distribution (Fig. 7b versus Fig. 7a). We therefore 
adopt these modifications for our baseline simulation 
(denoted as case specification B in Table 2; aside we 
note that retaining case specification A throughout 
would not alter our conclusions). Finally, considering 
this baseline simulation, we note that the simulated 
reflectivity is dominated by particles of 500-1000 gm in 
maximum dimension (Fig. 9a). 


d. Additional IN required to match observations 

Ice properties are quite uniform vertically in the base- 
line simulation (Fig. 10), as commonly observed (e.g., 
McFarquhar et al. 2007, 2011), including the total 
concentration of ice crystals (V,-, to which particles smaller 
than 200 gm contribute little (cf. Fig. 9a). In addition, 
N t below cloud, which is representative of the whole 
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Maximum Dimension (pm) 



Fig. 9. Simulated normalized contribution of ice to total 
dZIdlogD for the corresponding mean 


boundary layer, is about 200 times smaller than the 
overlying A^n of 1.7 L 1 (Table 3). To better match the 
mean observed ice crystal size distribution, we find that 
we need to initialize /V [N to a value 30 times greater than 
1.7 L 1 (IN X 30 in Fig. 7c and Tables 2 and 3). This also 
improves agreement with radar reflectivity (Fig. 8e), 
although the simulated range of radar reflectivity re- 
mains somewhat narrow and mean Doppler velocities 
somewhat slow (Figs. 8e and 8f), as discussed further 
below. When A IN is thus increased, ice crystal number 
concentration increases roughly linearly at all sizes, such 
that the normalized contributions of each particle size 
to number concentration and reflectivity remain nearly 
constant (Fig. 9b versus Fig. 9a). Despite the greater 
than tenfold increase in A,-, LWP develops only a modest 
downward trend and droplet concentration is negligibly 
impacted (Fig. 6, long dashed lines). 

The somewhat worsened agreement of simulated 
Doppler velocities with measurements prompts consid- 
eration of sensitivity to ice crystal habit, fall speed, and 
growth rate. Aggregates with plates fall faster than the 
radiating assemblages of plates assumed thus far (Fig. 11) 
and these two crystal types may be difficult to distinguish 
in some CPI images (see Fig. 4). Singular hexagonal 
plates fall slower and are relatively common in CPI 
images, but because their fall speeds are very similar to 
radiating plates over D of 100-400 /im, we use plates 



1 00 1 000 

Maximum Dimension (pm) 



number concentration dNIdlogD and radar reflectivity 
ice size distributions shown in Fig. 7. 


with sectorlike branches (sectored plates, also seen in 
CPI images) for a second sensitivity test. Assuming ag- 
gregates, simulated ice crystal size distributions shift to 
smaller sizes (Fig. 7d), radar reflectivity is correspond- 
ingly underestimated by about 10 dBZ (Fig. 8g), and 
mean Doppler velocities increase (Fig. 8h). By contrast, 
assuming sectored plates has a more modest, opposite 
effect (Figs. 7e and 8i,j). Size-resolved contributions to 
radar reflectivity shift accordingly (Figs. 9c,d). In gen- 
eral, given faster-falling crystals, more IN aloft would be 
required to match observed ice size distributions and 
radar reflectivities. 

We note that changes in the mode of IN activation 
have a lesser impact on our results than the foregoing 
changes in assumed ice habit (see Table 3), despite dif- 
ferences in the nucleated ice crystal size (e.g., D = 2 fim 
assumed for nucleated deposition IN versus preferen- 
tially large droplet size for nucleated immersion IN); 
when IN become available under cloud-top conditions, 
they are efficiently consumed in any mode except con- 
tact, and they grow rapidly to D > 100 gm regardless 
of initial size (cf. Fig. 9). Reducing the ratio of capaci- 
tance to maximum dimension to 0.35-0.25 for all crystal 
sizes (see section 3b) also has a lesser impact on size 
distribution shape (Fig. 7f versus Fig. 7c). However, the 
treatment of vapor growth rates for the diversity of ra- 
diating plates and other ice particle shapes seen in CPI 
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Fig. 10. Simulated profiles of domain-average ice nucleus number concentration N m , ice mass mixing ratio g,-, total ice crystal number 
concentration JV,-, and relative humidity over ice (RH,-) averaged over last 2 h of simulation time (corresponding to 2200-2400 UTC). 
Simulations listed in Table 2: baseline (solid lines), IN X 30 (dotted lines), aggregates (short dashed lines), plates (dash-dotted lines), 
modified diagnostic IN (dash-triple-dotted lines), and GCSS submission (long dashed lines). 


images is uncertain and should be considered further in 
future work. 

The sensitivity of results to cloud-top entrainment rate 
w e should also be considered. In simulations, w e (—0.1 
cm s , see Table 3) is computed as the rate of change of 

boundary layer depth H ( 0.3 cm s see Fig. 6) plus 

the large-scale subsidence rate at cloud top (—0.4 cm s 
see appendix), which is poorly constrained by reanalysis 
fields. However, in order to account for observed ice at 
quasi-equilibrium with above-cloud Nm of 1.7 L , w e 


would need to increase by a factor of 30 from about 0.1 to 
3 cm s '. It is difficult to reproduce relatively steady H 
under low-LWP conditions with such a large w e . For in- 
stance, if large-scale subsidence rate is increased by a 
factor of 2 and q v advective convergence increased suffi- 
ciently to maintain quasi-equilibrium LWP, then H de- 
creases by about 80 m over 4 h (not shown, compared with 
about 40 m in the baseline simulation and about 
70 m estimated from radar measurements), w e is reduced 
by about 20%, and more IN aloft would again be required 


Table 3. Simulation results: ice nucleus number concentration above the boundary layer N 1n , cloud-top entrainment rate w e , mean 
number-weighted ice crystal fall speed at the surface Vf, ice crystal number concentration predicted by Eq. (3) (Nmwjvf), ice crystal 
concentration JV,- and ice mass mixing ratio g, averaged below 280 m (representative of boundary layer values), and NJN in . All values 
averaged over the last 2 h of simulation time (corresponding to 2200-2400 UTC) except Nm, which is a model input (see Table 2). 


Simulation 

Nin (L *) 

w e (cm s x ) 

Vf (cm s *) 

Nmwjvf (L ') 

Ni (L- 1 ) 

Qi (mg kg ') 

Ni/N m (-) 

Steady-sate prognostic IN 

1.7 

0.17 

30.0 

0.0096 

0.0088 

0.025 

0.0052 

Baseline 

1.7 

0.13 

31.0 

0.0071 

0.0074 

0.021 

0.0043 

IN X 30 

51.0 

0.11 

30.0 

0.18 

0.29 

0.81 

0.0057 

Deposition IN only 

51.0 

0.11 

31.0 

0.18 

0.28 

0.77 

0.0055 

Condensation IN only 

51.0 

0.11 

30.0 

0.19 

0.26 

0.67 

0.0051 

Immersion IN only 

51.0 

0.12 

31.0 

0.20 

0.29 

0.81 

0.0057 

Decreased capacitance 

51.0 

0.12 

26.0 

0.24 

0.35 

0.55 

0.0069 

Aggregates 

51.0 

0.12 

38.0 

0.16 

0.23 

0.28 

0.0043 

Plates 

51.0 

0.12 

25.0 

0.24 

0.33 

1.2 

0.0065 

Modified diagnostic IN 

0.29* 

0.12 

27.0 

— 

0.32 

0.72 

— 

GCSS submission 

1.7* 

0.29 

32.0 

— 

1.8 

7.4 

— 


* Small differences between diagnostic N m and below-cloud mean V, attributable to boundary layer mixing conserving mixing ratio 
rather than concentration. 
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Fig. 11. Simulation ice particle fall speeds versus maximum di- 
mension calculated at the surface per Bohm (1989, 1999) for ra- 
diating plates (baseline, solid curves), aggregates (dashed curves), 
and plates with sectorlike branches (dash-dotted curves). Using 
the same ice crystal properties (see Table 1), fall speeds calculated 
per Heymsfield and Westbrook (2010) are shown for comparison 
(dotted curves). 

to match observations. We therefore believe that a factor 
of 30 increase in IN concentration likely errs on the low 
side required to explain the average in situ and remote 
sensing measurements within this modeling framework. 
Aside we note that even if local boundary layer depth 
were stationary, mesoscale gradients in boundary layer 
depth could exist that would not be captured with peri- 
odic boundary conditions (e.g., Avramov and Harrington 
2010). Lacking reliable observations of regional bound- 
ary layer depth gradients, this possibility is not pursued 
here. 

e. Entrainment limitations on rapidly consumed IN 

In all simulations with prognostic IN, equilibrium A* 
is two orders of magnitude smaller than A IN overlying 
the boundary layer (see Table 3). To understand the 
processes controlling A, -/Ain in these simulations, it is 
useful to consider a simple mixed-layer model for A, in 
the cloud-topped boundary layer, using the framework 
developed by Lilly (1968). As described above, ice 
properties are quite uniform vertically, entrained IN 
are rapidly activated, no other ice formation process 
is active, and the sole fate of all ice crystals is sedimen- 
tation to the surface. For a boundary layer of depth H 
entraining overlying air at a rate w e , ice crystals are 
therefore added at a rate wyA IN /// and sedimented at a 
rate VfNjIH , where ty is the number-weighted ice crystal 
fall speed at the surface. Cloud-top entrainment of ice- 
free air also dilutes A,- at a rate w e NjlH. Neglecting 
large-scale horizontal advective tendencies and the 
vertical dependence of air density, the mixed-layer A, 
budget can then be expressed as 


= w e N m - ( v f + w e) N i- (!) 

For the simulated conditions, w e iy (see Table 3), and 
Eq. (1) can be simplified to 

dN- 

= ^ 1N - yv (2) 

Dividing the ice crystal reservoir / /A, by its sink iy/V,- 
gives an e-folding time scale Hlvfoi about 20-30 min on 
which Nj relaxes toward its steady-state value 

N i = N W W e lv f ( 3 ) 

Table 3 shows A IN , w e , and iy averaged over hours 2-4 
of simulation time (cf. Fig. 6), the solution to Eq. (3), A,- 
averaged over hours 2-4 below cloud (representative 
of mean boundary layer values; cf. Fig. 10) and the ratio 
A//Vin- Equation (3) reproduces A,- to within 10% at the 
lower A, values and to within 30%-40% at the higher A,- 
values, in all cases capturing the two orders of magni- 
tude difference between A, : and A IN . Thus, A,/A IN -C 1 
since the supply of IN to the boundary layer is limited by 
a cloud-top entrainment rate that is much smaller than 
the number-weighted ice crystal fall speed. 

Quasi-equilibrium A, can be reached in simulations 
because the 4-h simulation time is about 10 times greater 
than the N, relaxation time. However, the divergence of 
agreement between Eq. (3) and simulated A,- at the 
higher A, values could be attributable to departure from 
quasi-equilibrium as desiccation increases (cf. Fig. 6) 
and cloud-top entrainment rate is reduced (cf. Table 3), 
consistent with ice loss rates that exceed supply rates in 
those simulations (Fig. 12). Equation (3) nonetheless 
explains the vast discrepancy between A,- and A IN . In 
addition, when Ain increases, the ice crystal size distri- 
bution shape remains relatively unaffected as it is shifted 
upward to greater A, (equivalent to multiplying by a 
size-independent factor), as noted above. Therefore tyis 
relatively constant, which is associated with a charac- 
teristic size distribution of ice in the boundary layer and 
a linear scaling of boundary layer ice mass mixing ratio 
< 7 , with A,- (Fig. 12). As shown above, the characteristic 
size distribution depends strongly on habit (e.g., as- 
suming aggregates rather than radiating plates reduces 
< 7 , by more than half), consistent with results from other 
case studies (e.g., Morrison and Pinto 2006; Avramov 
and Harrington 2010). 

f Modified diagnostic IN 

It is worthwhile to briefly compare our results with the 
GCSS SHEBA model intercomparison study (Morrison 
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Fig. 12. Simulation results with prognostic IN (diamonds) and diagnostic IN (triangles) averaged over the last 
2 h of simulation time (corresponding to 2200-2400 UTC). (left) Prognostic IN only: supply rate of IN to the 
boundary layer ./V 1 n H',, vs loss rate of ice crystals to the surface /V,iy, where dashed line indicates 1:1. (right) Domain- 
mean Nj vs q t below cloud base, where dashed line indicates linear relation in simulations with prognostic IN and 
radiating plates, and dotted lines indicate linear relations expected with prognostic IN and aggregates (lower dotted 
line) or sectored plates (upper dotted line). Values taken from Table 3 for the following simulations: steady-state 
prognostic IN, baseline, IN X 30, aggregates, plates, modified diagnostic IN, and GCSS submission. 


et al. 2011), where DHARMA ice properties were based 
on radiating plates as in most simulations here. A near- 
equilibrium LWP was achieved with 1.7 L _1 IN treated 
diagnostically in the DHARMA baseline submission to 
that study (Morrison et al. 2011, their Fig. 4), in contrast 
to the rapid loss of LWP found in this study (see Fig. 6, 
solid lines). This is principally because the specified 
horizontal advective moisture convergence, which gen- 
erally cannot be adequately constrained by reanalysis 
fields, was larger in the GCSS SHEBA case (see Fig. A4) 
and was therefore able to balance a higher rate of des- 
iccation associated with TV, of 1.7 L 1 . In simulations of 
mixed-phase Arctic clouds observed during the Beaufort 
Arctic Storms Experiment, Jiang et al. (2000) demon- 
strated how an observed quasi-equilibrium LWP can be 
achieved over a wide range of possible TV,- when offsetting 
changes in advective tendencies are made. The very large 
uncertainty in observations of both TV, and advective 
tendencies therefore introduces a large corresponding 
uncertainty in modeling case studies. 

However, in our baseline DHARMA submission to 
the GCSS SHEBA intercomparison study, we found 
that radar reflectivities during 1200-2400 UTC exceeded 
MMCR measurements by greater than 10 dBZ (Fig. 8m), 
consistent with the possibility that ice crystal number 
concentrations were too high (Fig. 7h). That the mean 
Doppler velocity distribution nonetheless appeared quite 
similar to measurements (Fig. 8n) suggested that simu- 
lated ice properties might be reasonable. Aside we note 
that LWP fell roughly fourfold over 1200-2400 UTC 


(Morrison et al. 2011, their Fig. 4), although radar re- 
flectivity and mean Doppler velocity distributions appear 
roughly similar during 1200-2400 and 2200-2400 UTC 
(see Figs. 8m, n). In this study, using the adjusted case 
specification, which achieves quasi-equilibrium LWP un- 
der ice-free conditions, we are able to simultaneously 
reproduce sustained LWP, ice crystal number size dis- 
tribution features, radar reflectivities, and mean Doppler 
velocities. However, this can only be done with prognostic 
IN, which always produces TV,- -C TV IN , and it also requires 
overlying TV IN to be elevated by a factor of about 30. Then 
simulated TV, reaches about 0.3 L , similar to the GCSS 
SHEBA study sensitivity test with diagnostic TV IN of about 
0.2 L 1 in which DHARMA and other models predict 
increased LWP. 

Finally, we find here that a modified diagnostic TV IN 
fixed at 0.29 L 1 rather than 1.7 L _1 (see Fig. 10) can 
also reproduce measurements quite well. Crystals smaller 
than 200 gm are enhanced below cloud compared with 
prognostic IN, but the size distribution of larger ice is 
minimally affected (Fig. 7g), leading to little change in 
radar reflectivities and mean Doppler velocities (Figs. 8k,l). 
Aside we note that had ice instead been treated as ag- 
gregates in our baseline submission to the GCSS SHEBA 
intercomparison study, median reflectivity would have 
dropped about 5 dBZ (not shown) versus dropping about 
10 dBZ with prognostic IN in this study; the sensitivity to ice 
habit using diagnostic IN is less than that using prognostic 
IN because the impact of habit on ice loss rate is flexibly 
compensated by an unlimited source of new ice crystals. 


January 2012 


FRIDLIND ET AL. 


381 



800 

E, 

600 

5: 

400 

o> 


'(D 

X 

200 


(a) Observed 



ol 

22.0 


22.5 


23.0 23.5 

Time (UTC) 


800 

1 600 
2 400 

O) 

X 200 
0 


dBZ 

0. 


800 

-10. 

E 

600 

-20. 

£ 

400 

-30. 

O) 


-40. 

'<D 

X 

200 

-50. 


0 


(b) Observed 


sHIF'i'T'* i 


m s 
2.0 
1.5 
1.0 
0.5 
0.0 
-0.5 


24.0 


22.0 


22.5 


23.0 

Time (UTC) 


23.5 


24.0 


(c) Baseline 

- 

"vnufliNir: 

nun 



dBZ 

0. 

-10. "E 
- 20 . ~ 


800 

600 


400 

-30. - 
-40. X 200 
-50. 0 


(d) Baseline 


riramm™ 


m s’ 1 
2.0 
1.5 
1.0 
0.5 
0.0 
-0.5 


100 


200 300 400 

Model Gridbox 


500 600 


100 


200 300 400 

Model Gridbox 


500 600 


800 

1 600 

£ 400 
o> 

X 200 
0 


(e) IN x 30 

- 

fimflfllitlMlN 




dBZ 

0. 

- 10 . 

- 20 . 

-30. 

-40. 

-50. 


800 
g 600 
£ 400 

O) 

X 200 
0 



m s' 1 
2.0 
1.5 
1.0 
0.5 
0.0 
-0.5 


0 100 200 300 400 500 600 

Model Gridbox 


0 100 200 300 400 500 600 

Model Gridbox 


Fig. 13. Observed and simulated (left) radar reflectivity and (right) mean Doppler velocity. Observed 35-GHz 
reflectivity and Doppler velocity measured by the MMCR during (a),(b) 2200-2400 UTC 7 May 1998. Simula- 
tions listed in Table 2: (c),(d) baseline and (e),(f) IN X 30. Simulation results calculated at 3 h (corresponding to 
2300 UTC). 


We note that solving Eq. (3) for ;V IN required to sup- 
port Nj of 1.7 L _1 using w e 0.3 cm s 1 and ty 
30 cm s 1 (see GCSS submission in Table 3) gives an /V| N 
of about 200 L ', which is high compared with typically 
measured conditions (e.g., DeMott et al. 2010). Overall, 
based on our model results compared with forward- 
simulated radar variables and in situ ice crystal size dis- 
tributions, we hypothesize that actual TV,- were sustained 
closer to 0.17 than 1.7 L . However, the discrepancy 
between /V IN observed and ;V IN required to reproduce 
observed ice properties in our adjusted case study here 
indicates that substantial problems remain in either the 
model, the case study formulation, and/or the observa- 
tional dataset. These results are rather similar to past 
findings in the 10 October M-PACE case study (e.g., 
Fridlind et al. 2007; Fan et al. 2009), but contrast with 
relatively greater success matching simultaneously ob- 
served N t and 7V rN in the 8 April ISDAC case study 
(Avramov et al. 2011). 

g. Horizontal heterogeneity of ice 

As shown above, when simulations approximately 
reproduce observed ice crystal size distributions, they 
also roughly reproduce observed radar reflectivities. 
However, radar observations indicate a range in reflec- 
tivity over 2200-2400 UTC that is notably greater than 
simulated (e.g., Fig. 8e). Furthermore, periods of low 
reflectivity were of extended duration (Fig. 13a), as 


during 2300-2350 UTC, which at about 30 min duration 
at cloud-level horizontal wind speeds corresponds to a 
horizontal distance of about 8 km that is about 10 times the 
boundary layer depth. Periods of similar duration were 
characterized by higher radar reflectivities. Observations 
during 1200-2200 UTC indicate that such variability was 
commonplace in this cloud deck (cf. Morrison et al. 2011). 

Using the visualization method described by van 
Diedenhoven et al. (2009), Figs. 13c and 13e illustrate 
how our simulations fail to reproduce observed vari- 
ability in reflectivity. Increasing domain size to 12.8 X 
12.8 km 2 produces indistinguishable results (not shown), 
consistent with a weak feedback of the nonsublimating ice- 
phase precipitation on convective dynamics, which con- 
trasts with a strong feedback of evaporating liquid-phase 
precipitation (cf. Feingold et al. 2010). But we cannot 
rule out that a much larger domain size or longer-duration 
case study would produce other results. Such pronounced 
alternating reflectivity features on horizontal scales many 
times greater than the boundary layer depth were not 
present in the 10 October M-PACE or 8 April ISDAC 
cases, where the observed variability of radar reflectivity 
was reliably reproduced by various simulations (van 
Diedenhoven et al. 2009; Avramov et al. 2011). Here 
periods of lower reflectivity tend to resemble the baseline 
simulation, whereas periods of higher reflectivity re- 
semble the simulation with enhanced IN concentration. It 
is uncertain what modifications to the model setup used 
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here would be needed to reproduce the observed degree 
of horizontal variability in cloud ice. 

5. Conclusions and implications 

We adjusted the GCSS SHEBA case study for mixed- 
phase boundary layer clouds observed during 1200- 
2400 UTC on 7 May 1998 (Morrison et al. 2011) in order 
to more closely match conditions during the briefer 
2000-2400 UTC time span when airborne ice particle 
size distribution measurements were obtained. Our prin- 
cipal objective is to determine whether simulations can 
reproduce all available measurements when using the 
mean ice nucleus (IN) number concentration A IN mea- 
sured above cloud. Since in situ measurements of ice 
crystal total number concentration N t were unreliable, 
we compare simulation results with 1) in situ measure- 
ments of the size distribution of ice with maximum di- 
mension larger than 200 /l m and 2) ground-based remote 
sensing measurements of cloud radar reflectivity and 
mean Doppler velocity. Results can be briefly summa- 
rized as follows. 

1) When A IN is initialized to the observed mean, treat- 
ing IN prognostically (accounting for consumption 
when activated) gives dramatically different results 
than treating IN diagnostically (neglecting consump- 
tion by definition), which is not a new finding (e.g., 
Harrington and Olsson 2001; Rasmussen et al. 2002; 
Morrison et al. 2005). Consumption depletes rapidly 
activated IN from a well-mixed boundary layer within 
minutes. This large difference in model results has 
important implications for interpreting simulations 
that treat IN diagnostically (e.g., Jiang et al. 2000). 
Namely, diagnostic A IN should be interpreted as in- 
cloud A,-, which may differ substantially from /\j N in 
cloud-free air that is entrained. 

2) When treating IN prognostically, simulated consump- 
tion proceeds rapidly in all nucleation modes except 
contact, which proceeds too slowly to be a significant 
source of ice crystals if IN are assumed to be 0.5 /im in 
diameter, the mean effective dimension measured by 
the CFDC during SHEBA (Rogers et al. 2001). It has 
been argued that contact nucleation could play an 
important role under SHEBA conditions (Morrison 
et al. 2005), but available measurements are insuffi- 
cient to constrain actual contact nucleation rates. 
Results are insensitive to the whether IN are alterna- 
tively activated in the deposition, condensation, or 
immersion modes. Here we have neglected possible 
effects of nucleation mode on ice habit (Bailey and 
Hallett 2002; Bacon et al. 2003), which could conceiv- 
ably be important since results are sensitive to habit. 


3) If rapidly activated IN are the principal source of new 
ice crystals, as commonly assumed (e.g., Fan et al. 
2009), we find that a factor of about 30 greater A IN 
than observed is required to reproduce observed ice 
crystal size distributions and cloud radar reflectivities 
when accounting for IN consumption. Although radar 
reflectivities are weighted toward larger particles than 
ice number size distributions, both exhibit peaks in the 
200-1000-q.m size range spanned by a single ice mode 
(cf. Figs. 7 and 9). Thus, measured A IN appear in- 
sufficient to explain observed ice in this case study. 
It is unknown to what degree the factor of about 
30 discrepancy found here can be attributed to 
observational uncertainties or modeling shortcomings. 
For instance, the CFDC is not designed to measure IN 
larger than about 2 /mi in diameter and may un- 
dercount IN active in the contact mode (e.g., Rogers 
et al. 2001; McFarquhar et al. 2011). In the 10 October 
M-PACE case, similar results using two independent 
models led to the speculative consideration of novel 
ice formation mechanisms unconstrained by CFDC 
measurements of A IN (e.g., Fridlind et al. 2007; Fan 
et al. 2009), but an 8 April ISDAC case study shows 
less discrepancy, which is on the order of experimental 
uncertainty (Avramov et al. 2011). Although blowing 
snow is not generally lifted at the low horizontal wind 
speeds observed in this case (e.g., Walden et al. 2003), 
it cannot be ruled out as a conceivable ice crystal source 
over pack ice. 

4) When IN are rapidly consumed, /V, is always more 

than two orders of magnitude smaller than overlying 
A IN . Under the microphysically simple conditions 
of this case (apparently negligible ice aggregation, 
sublimation, and multiplication), in the equilibrium 
state of a simple mixed-layer model for A,- [Eq. (3)], 
A,-/A| N equals the entrainment rate w e divided by 
the number-weighted ice fall speed at the surface ly. 
Here uy/iy -C 1 since vty < 1 cm s _1 and iy > 
10 cm s . Conditions where A,-/A IN 1 contrast 
with conditions where the IN supply rate is not limited 
by entrainment (e.g., in a wave cloud; Eidhammer 
et al. 2010), with implications for interpreting regional 
measurements. For instance, Prenni et al. (2009) point 
to observations of A, A IN in the Arctic as evidence 
that 1) observed A in are adequate to explain observed 
ice and 2) secondary ice sources are not important. 
But here A,- A IN would be evidence that secondary 
ice sources must be important (otherwise A ; <C A IN ). 
Finally, to the extent that A;/ Ain ny/iy < 1, en- 

trainment rate limitations on the IN supply rate serve as 
a buffer on cloud system sensitivity to increasing over- 
lying A in in the sense outlined by Stevens and Feingold 
(2009). Aggregation could decrease sensitivity, whereas 
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multiplication (e.g., Yano and Phillips 2011) could 
increase it. We note that blowing snow, seeding from 
aloft, or any other ice crystal formation processes not 
related to entrained IN or existing ice would in- 
troduce independent source terms in Eq. (1). 

5) Simulations fail to reproduce the observed horizon- 
tal heterogeneity of radar reflectivity even when 
domain size is increased. Pronounced alternating 
increases and decreases of reflectivity on horizontal 
scales of about 10 times the boundary layer depth 
distinguish this case from the 10 October M-PACE 
and 8 April ISDAC cases with mixed-phase cloud 
layers (see section 4g). The contributing dynamical 
and microphysical causes are unknown, and it is 
uncertain what modifications to the model setup 
used here would be needed to reproduce the observed 
heterogeneity. 

6) Simulation results are sensitive to assumed ice prop- 
erties not adequately constrained by measurements, 
consistent with previous work (e.g., Morrison and 
Pinto 2006; Avramov and Harrington 2010). The 
irregular habits that exist in mixed-phase clouds pres- 
ent a challenge to models (e.g., Bailey and Hallett 
2002). For this study, the most appropriate observa- 
tional constraints would have been direct single- 
particle field measurements of ice crystal mass, 
maximum dimension, projected area, aspect ratio, 
and terminal fall speed, suitable to identify both 
mean properties and their spread. Such measure- 
ments could be made simultaneously at ground level 
(e.g., Kajikawa 1972), perhaps in part by instruments 
that could be deployed unattended (e.g., Newman 
et al. 2009; Barthazy et al. 2004). 
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Fig. Al. Observed and simulated domain-mean LWP, surface 
downwelling shortwave and longwave radiative fluxes (SW down and 
LW down ), surface sensible and latent heat fluxes (SHF and LHF), 
and boundary layer depth H (defined by mean elevation where 
liquid water potential temperature is 258 K; not recorded in ob- 
servations). Simulations: bulk warm microphysics with the GCSS 
model intercomparison specification (solid lines), with surface 
fluxes predicted using similarity theory (dotted lines), and with 
adjusted initial conditions and large-scale forcings (dashed lines). 
Observed range (shaded) is during 2200-2400 UTC 7 May 1998 
with estimated uncertainty (cf. Persson et al. 2002). 

NCEP reanalysis data provided by the NOAA/OAR/ 
ESRL PSD, Boulder, Colorado, from their website at 
http://www.esrl.noaa.gov/psd/. 

APPENDIX 

Case Study Development 

Here our objective is to make several adjustments 
to the 1200-2400 UTC 7 May 1998 GCSS SHEBA case 
study (Morrison et al. 2011) in order to achieve 
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Fig. A2. Observed and simulated profiles of temperature T, water vapor mixing ratio q v , potential temperature 6, horizontal wind speed 
(WS), and relative humidity (RH). Observations at 1800 and 2400 UTC (plus and asterisk symbols, respectively). Simulations using bulk 
warm microphysics with the GCSS model intercomparison specification (solid lines at 0 h, initial condition) and with adjusted initial 
conditions and large-scale forcings (dashed and dash-dotted lines at 0 and 4 h), and in baseline simulations using mixed-phase bin 
microphysics (dotted line at 4 h, initial condition same as adjusted case) and with IN X 30 (dash-triple-dotted lines at 4 h, initial condition 
same as adjusted case). Simulation times of 0 and 4 h correspond to 2000 and 2400 UTC. 


reasonably close simultaneous agreement with the 
following observed conditions specifically during 2200- 
2400 UTC: liquid water path (LWP), surface upwelling 
and downwelling radiative fluxes, surface skin tem- 
perature, 10-m tower measurements of temperature 
and water vapor, surface turbulent heat fluxes, and 
observed profiles of temperature, water vapor, poten- 
tial temperature, and wind speed. The GCSS model 
intercomparison specification for 12-h simulations was 
based on a combination of observations, reanalysis 
fields (European Centre for Medium-Range Weather 
Forecasts), and model results. Here we shorten the 
simulation time to 4 h, allowing 2 h for model spinup 
before comparison of simulated conditions during 
hours 2-4 with observations from 2200 to 2400 UTC. 

The effects of ice can be neglected during case study 
development because desiccation remains relatively 
weak when ice crystal number size distribution features 
match the available observations (shown in section 4). 
We therefore save computational time by using an effi- 
cient bulk warm microphysics scheme that consists of 
condensational adjustment with slow sedimentation of 
cloud droplets following Ackerman et al. (2009). Droplet 
number concentration N,/ is fixed at 215 cm , consistent 
with observations. At the very low ratio of observed LWP 
to Nd here (—0.02 g in 2 cm 3 ), gravitational collection can 
be neglected (cf. Comstock et al. 2004, their Fig. 10). A 
lognormal droplet size distribution with a geometric 
standard deviation of 1.3 is assumed for radiative transfer 
and sedimentation. 


We start with the initial and boundary conditions and 
large-scale forcings from the GCSS SHEBA case. LWP 
is initially 20 g m 2 . consistent with an average of obser- 
vations over 1200-2400 UTC but greater than observed 
during 2200-2400 UTC (Fig. Al, solid lines). At 4 h, 
simulated LWP is roughly 5 times greater than observed, 
resulting in underprediction and overprediction of surface 
downwelling shortwave and longwave radiative fluxes, 
respectively. Predicted LWP is not sensitive to replacing 
the fixed surface latent and sensible heat fluxes specified 
in the model intercomparison with interactive fluxes 
predicted at grid scale using similarity theory (Fig. Al, 
dotted lines). 

We first make several adjustments to reduce initial 
LWP and simultaneously improve consistency with 2200- 
2400 UTC observations. The initial temperature T profile 
is made uniformly colder by 0.5 K, the water vapor mixing 
ratio in the boundary layer limited to q v < 0.829 g kg~\ 
and all profiles shifted downward by 10 m. These changes 
bring initial conditions closer to the 1800 UTC sounding 
(Fig. A2). Aside we note that reported q„ and T profiles 
correspond to an LWP that is far greater than retrieved 
from observations, presumably owing to measurement 
bias. We accept T as the more reliably measured pa- 
rameter and use reported LWP to constrain initial q v . 
Surface skin temperature is increased by 1 K (Fig. A3). 

With adjustments to initial and boundary conditions 
in place, we turn next to large-scale forcing terms. To 
emulate the observed evolution of the T profile, we in- 
crease the potential temperature horizontal advective 
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Fig. A3. Observed and simulated time series of LWP, upwelling and downwelling shortwave and longwave ra- 
diative fluxes (SW down , SW up , LW down , LW llp ), surface skin temperature, 10-m wind speed (WS), 10-m air temper- 
ature and water vapor mixing ratio q v , surface SHF and LHF, and 10-m LHF. Hourly observations calculated from 
tower and surface measurements [asterisks with uncertainty range (cf. Persson et al. 2002), see section 2], Simulations 
use bulk microphysics with adjusted initial conditions and large-scale forcings (dashed lines), and mixed-phase bin 
microphysics in the baseline case (dash-dotted lines) and with IN X 30 (dotted lines). 


tendency to a vertically uniform value of 2 K day -1 . This 
is larger at most elevations than in the 12-h GCSS speci- 
fication and larger than indicated from analysis of NCEP 
fields (Fig. A4), but we consider agreement with 1800 and 
2400 UTC soundings a better constraint. Last, we adjust 
the q v horizontal advective tendency to a vertically uni- 
form value of 0.02 g kg -1 day -1 in order to achieve LWP 
within the observed range and at quasi-equilibrium, which 
we define for a parameter by requiring its e-folding time 
(computed from a 1-h running mean of domain averages 
reported every minute) to continuously exceed 10 h during 
simulation hours 3 and 4. The adjusted moisture hori- 
zontal advective tendency is smaller than most NCEP 
values, but is constrained relative to other forcings if LWP 
is to maintain quasi-equilibrium (see also section 4f). 


A simulation with all foregoing changes achieves rel- 
atively close agreement with LWP, upwelling and 
downwelling surface radiative fluxes, surface skin T, 10-m 
tower measurements of T and <y„, surface sensible heat 
flux, and observed profiles of T, q v , potential temper- 
ature, wind speed, and relative humidity (see Figs. A2 
and A3). A notable exception is disagreement with 
measured surface latent heat flux, although agreement 
with bulk calculations at 10 m is very good; the cause 
for persistent disagreement between bulk calculations 
and measured latent heal fluxes at the surface is un- 
known (Persson et al. 2002). Finally, a net effect of all 
adjustments is a reduction in cloud-top entrainment 
such that boundary layer depth falls by about 30 m over 
4 h (see Fig. Al), which improves consistency with 
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Fig. A4. Derived profiles of large-scale vertical wind w L s and advective tendencies of water vapor and potential 
temperature specified in the GCSS intercomparison (solid lines) and derived from NCEP reanalysis fields for 7 May 
1998 (other line types, see legend). 


MMCR observations and soundings somewhat (cf. 
Fig. A2). 
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